1 Intro to Spatial Data

In the workshop, we learned about:

  • What is Spatial Data?
  • What is the sf framework for R?

To delve in further, let’s see some spatial data in action.

We’ll work with the sf library first.

1.1 Load Spatial Data

First load in the shapefile. Remember, this type of data is actually comprised of multiple files. All need to be present in order to read correctly.

## Reading layer `geo_export_aae47441-adab-4aca-8cb0-2e0c0114096e' from data source 
##   `/Users/maryniakolak/code/Intro2RSpatialMed/data/geo_export_aae47441-adab-4aca-8cb0-2e0c0114096e.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 801 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.94025 ymin: 41.64429 xmax: -87.52366 ymax: 42.02392
## CRS:           4326

1.2 Non-Spatial & Spatial Views

Always inspect data when loading in. First we look at a non-spatial view.

## Simple feature collection with 6 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.68822 ymin: 41.72902 xmax: -87.62394 ymax: 41.87455
## CRS:           4326
##   commarea commarea_n countyfp10     geoid10 name10        namelsad10
## 1       44         44        031 17031842400   8424 Census Tract 8424
## 2       59         59        031 17031840300   8403 Census Tract 8403
## 3       34         34        031 17031841100   8411 Census Tract 8411
## 4       31         31        031 17031841200   8412 Census Tract 8412
## 5       32         32        031 17031839000   8390 Census Tract 8390
## 6       28         28        031 17031838200   8382 Census Tract 8382
##   notes statefp10 tractce10                       geometry
## 1  <NA>        17    842400 POLYGON ((-87.62405 41.7302...
## 2  <NA>        17    840300 POLYGON ((-87.68608 41.8229...
## 3  <NA>        17    841100 POLYGON ((-87.62935 41.8528...
## 4  <NA>        17    841200 POLYGON ((-87.68813 41.8556...
## 5  <NA>        17    839000 POLYGON ((-87.63312 41.8744...
## 6  <NA>        17    838200 POLYGON ((-87.66782 41.8741...

Note the last column – this is a spatially enabled column. The data is no longer a ‘shapefile’ but an `sf’ object, comprised of polygons.

We can use a baseR function to view the spatial dimension. The sf framework enables previews of each attribute in our spatial file.

1.3 Spatial Data Structure

Check out the data structure of this file… What object is it?

## Classes 'sf' and 'data.frame':   801 obs. of  10 variables:
##  $ commarea  : Factor w/ 77 levels "1","10","11",..: 39 55 28 25 26 21 62 49 74 75 ...
##  $ commarea_n: num  44 59 34 31 32 28 65 53 76 77 ...
##  $ countyfp10: Factor w/ 1 level "031": 1 1 1 1 1 1 1 1 1 1 ...
##  $ geoid10   : Factor w/ 801 levels "17031010100",..: 785 767 772 773 756 751 584 513 684 34 ...
##  $ name10    : Factor w/ 801 levels "1001","1002",..: 782 764 769 770 753 748 545 443 663 266 ...
##  $ namelsad10: Factor w/ 801 levels "Census Tract 1001",..: 782 764 769 770 753 748 545 443 663 266 ...
##  $ notes     : Factor w/ 7 levels "Half in CA 64 (Midway Airport)",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ statefp10 : Factor w/ 1 level "17": 1 1 1 1 1 1 1 1 1 1 ...
##  $ tractce10 : Factor w/ 801 levels "010100","010201",..: 785 767 772 773 756 751 584 513 684 34 ...
##  $ geometry  :sfc_POLYGON of length 801; first list element: List of 1
##   ..$ : num [1:243, 1:2] -87.6 -87.6 -87.6 -87.6 -87.6 ...
##   ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA
##   ..- attr(*, "names")= chr  "commarea" "commarea_n" "countyfp10" "geoid10" ...

Check out the coordinate reference system. What is it? What are the units?

## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS84(DD)",
##     DATUM["WGS84",
##         SPHEROID["WGS84",6378137.0,298.257223563]],
##     PRIMEM["Greenwich",0.0],
##     UNIT["degree",0.017453292519943295],
##     AXIS["Geodetic longitude",EAST],
##     AXIS["Geodetic latitude",NORTH]]

1.4 Exploring Coordinate Reference Systems

Lets see how switching CRS changes our object. First we’ll try the Mollweide coordinate reference system that does a good job preserving area across the globe.

To transform our CRS, we use the st_transform function. To plot, we use baseR again but with some paremeter updates. Finally, we check out the CRS of our new object. What are the units? Any other details to note? Will this be appropriate for our spatial analysis?

## Coordinate Reference System:
##   User input: +proj=moll 
##   wkt:
## PROJCS["unnamed",
##     GEOGCS["WGS 84",
##         DATUM["unknown",
##             SPHEROID["WGS84",6378137,298.257223563]],
##         PRIMEM["Greenwich",0],
##         UNIT["degree",0.0174532925199433]],
##     PROJECTION["Mollweide"],
##     PARAMETER["central_meridian",0],
##     PARAMETER["false_easting",0],
##     PARAMETER["false_northing",0]]

Next, we’ll try the Winkel CRS, which is a compromise projection that facilitates minimal distortion for area, distance, and angles. We use the same approach, recyling the code with new inputs.

## Coordinate Reference System:
##   User input: EPSG:54019 
##   wkt:
## PROJCS["World_Winkel_II",
##     GEOGCS["GCS_WGS_1984",
##         DATUM["WGS_1984",
##             SPHEROID["WGS_84",6378137.0,298.257223563]],
##         PRIMEM["Greenwich",0.0],
##         UNIT["Degree",0.0174532925199433]],
##     PROJECTION["Winkel_II"],
##     PARAMETER["False_Easting",0.0],
##     PARAMETER["False_Northing",0.0],
##     PARAMETER["Central_Meridian",0.0],
##     PARAMETER["Standard_Parallel_1",50.45977625218981],
##     UNIT["Meter",1.0],
##     AUTHORITY["Esri","54019"]]

We could also try a totally different projection, to see how that changes our spatial object. Let’s use the “Old Hawaiian UTM Zone 4n” projection, with the EPSG identified from an online search. How does this fare?

Finally.. let’s choose a projection that is focused on Illinois, and uses distance as feet or meters, to make it a bit more accessible for our work. EPSG:3435 is a good fit:

## Coordinate Reference System:
##   User input: EPSG:3435 
##   wkt:
## PROJCS["NAD83 / Illinois East (ftUS)",
##     GEOGCS["NAD83",
##         DATUM["North_American_Datum_1983",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6269"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4269"]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",36.66666666666666],
##     PARAMETER["central_meridian",-88.33333333333333],
##     PARAMETER["scale_factor",0.999975],
##     PARAMETER["false_easting",984250.0000000002],
##     PARAMETER["false_northing",0],
##     UNIT["US survey foot",0.3048006096012192,
##         AUTHORITY["EPSG","9003"]],
##     AXIS["X",EAST],
##     AXIS["Y",NORTH],
##     AUTHORITY["EPSG","3435"]]

1.5 Refine Basic Map

Now we’ll switch to a more extensive cartographic mapping package, tmap. We approach mapping with one layer at a time. Always start with the object you want to map by calling it with the tm_shape function. Then, at least one descriptive/styling function follows. There are hundreds of variations and paramater specifications, so take your time in exploring tmap and the options.

Here we style the tracts with some semi-transparent borders.

Next we fill the tracts with a light gray, and adjust the color and transparency of borders. We also add a scale bar, positioning it to the left and having a thickness of 0.8 units, and turn off the frame.

Check out https://rdrr.io/cran/tmap/man/tm_polygons.html for more ideas!

1.6 Arrange multiple maps

Sometimes we want to look at multiple maps at once. Write your mapping function to a new variable, and then call that variable in order of desire using the tmap_arrange function. Hint: this is just one of many! ways to map multiples using tmap… see if you can uncover more in the documentation.

1.7 Interactive Mode

So far, we’ve been plotting static maps. We can also switch to an interactive map that uses a Leaflet widget by switching the tmap_mode() parameter specification from “plot” to “view.” It’s on “plot” as default.

## tmap mode set to interactive viewing

Map the same map as before, and check out the interaction!

The tracts are not transparent enough, so we update that here. You can also click the box on the left side to try out other basemaps. See if you can find out how to add a basemap to a static/plotted map, using tmap documentation…

We revert back to plot mode for now.

## tmap mode set to plotting

1.8 Overlay Zip Code Boundaries

How do census tract areas correspond to zip codes? While tracts better represent neighborhoods, often times we are stuck with zip code level scale in healh research. Here we’ll make a reference map to highlight tract distribution across each zip code.

First, we read in zip code boundaries. This data was downloaded directly from the City of Chicago Data Portal as a shapefile.

## Reading layer `geo_export_54bc15d8-5ef5-40e4-8f72-bb0c6dbac9a5' from data source 
##   `/Users/maryniakolak/code/Intro2RSpatialMed/data/geo_export_54bc15d8-5ef5-40e4-8f72-bb0c6dbac9a5.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 61 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
## CRS:           4326

Next, we layer the new shape in – on top of the tracts. We use a thicker border, and try out a new color. Experiment!